Spotify Data Analysis and Descriptive Statistics

1. Data Set Description For this project, our team decided to take a data-driven approach to evaluate music. We are using the Spotify Dataset from 1921 to 2020 that consists of over 160,000 different tracks. In order to stay consistent in the evaluation, all data was sourced from the Spotify Web API. The following set of data is combination of Primary, Numerical, Dummy(binary), and Categorical data. This allows us to explore different types of models and draw different insights. Here are all the variables that is included in the data set:

Primary - id: this an unique key comprised of numbers and characters that is assigned to each track generated by Spotify

Numerical - acousticness: range from 0(LOW) to 1(HIGH) - danceability: range from 0(LOW) to 1(HIGH) - energy: range from 0(LOW) to 1(HIGH) - duration_ms: majority range from 200,000 to 300,000 - instrumentalness: range from 0(LOW) to 1(HIGH) - valence: range from 0(LOW) to 1(HIGH) - popularity: range from 0(LOW) to 100(HIGH)* - tempo: majority range from 50(LOW) to 150(high) - liveness: range from 0(LOW) to 1(HIGH) - loudness majority range from -60 to 0 - speechiness: range from 0(LOW) to 1(HIGH) - year: range from 1921 to 2020

Dummy - mode: 0 represents minor and 1 represents major - explicit: 0 represents no explicit content and 1 represents explicit content

Categorical - key: this consists of all different music keys on onctave encoded from 0 to 11 (i.e. C = 0, C# = 1, etc…) - artists: the artist of the track - release_date: the date of release in yyyy-mm-dd format - name: the name of the track

*With our approach on evaluating different tracks, we decided to use popularity as our main dependant variable.

options(repr.matrix.max.rows=100, repr.matrix.max.cols=20)

Import all Libraries

# Libraries
options(warn=-1)
library(ggplot2)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(plotly)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

library(hrbrthemes)

NOTE: Either Arial Narrow or Roboto Condensed fonts are required to use these themes.
      Please use hrbrthemes::import_roboto_condensed() to install Roboto Condensed and
      if Arial Narrow is not on your system, please see https://bit.ly/arialnarrow

library(forecast)

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

library(xts)

Loading required package: zoo

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric


Attaching package: ‘xts’

The following objects are masked from ‘package:dplyr’:

    first, last

library(Metrics)

Attaching package: ‘Metrics’

The following object is masked from ‘package:forecast’:

    accuracy

library(psych)

Attaching package: ‘psych’

The following objects are masked from ‘package:ggplot2’:

    %+%, alpha

library(dygraphs)
library(GGally)

Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ tibble  3.0.3     ✓ purrr   0.3.4
✓ tidyr   1.1.2     ✓ stringr 1.4.0
✓ readr   1.4.0     ✓ forcats 0.5.0
── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x psych::%+%()     masks ggplot2::%+%()
x psych::alpha()   masks ggplot2::alpha()
x plotly::filter() masks dplyr::filter(), stats::filter()
x xts::first()     masks dplyr::first()
x dplyr::lag()     masks stats::lag()
x xts::last()      masks dplyr::last()

library(tidyquant)  

Loading required package: lubridate

Attaching package: ‘lubridate’

The following objects are masked from ‘package:base’:

    date, intersect, setdiff, union

Loading required package: PerformanceAnalytics

Attaching package: ‘PerformanceAnalytics’

The following object is masked from ‘package:graphics’:

    legend

Loading required package: quantmod
Loading required package: TTR
Version 0.4-0 included new data defaults. See ?getSymbols.
══ Need to Learn tidyquant? ═══════════════════════════════════════════════════════════════════════════════════════════════════════════════
Business Science offers a 1-hour course - Learning Lab #9: Performance Analysis & Portfolio Optimization with tidyquant!
</> Learn more at: https://university.business-science.io/p/learning-labs-pro </>

library(cranlogs)   
library(corrr)      
library(cowplot)  

Attaching package: ‘cowplot’

The following object is masked from ‘package:lubridate’:

    stamp

library(Metrics)

Initial data visualization

Plotting annual median values for each feature overall. Standardized such that each feature is on the same scale.

df <- read.csv("./data/data.csv")
q1s = data.frame(matrix(ncol=3, nrow=0))
medians = data.frame(matrix(ncol=3, nrow=0))
q3s = data.frame(matrix(ncol=3, nrow=0))
feature_names = c("acousticness", "danceability", "instrumentalness", "energy",
                  "duration_ms", "valence", "tempo", "liveness", "loudness", "speechiness")
for(i in 1921:2020){
  annual_val = subset(df, year==i)

  for(name in feature_names){
    vals = quantile(annual_val[,name])[2:4]
    q1s <- rbind(q1s, c(i, name, vals[1]))
    medians <- rbind(medians, c(i, name, vals[2]))
    q3s <- rbind(q3s, c(i, name, vals[3]))
  }
}

colnames(q1s) <-c("Year", "Feature", "Value")
colnames(medians) <-c("Year", "Feature", "Value")
colnames(q3s) <-c("Year", "Feature", "Value")

split_q1s <- split(q1s, q1s$Feature)
split_medians <- split(medians, medians$Feature)
split_q3s <- split(q3s, q3s$Feature)

feature_quantiles = array(0, c(100,3,length(feature_names)))

for(i in 1:length(feature_names)){
  for(j in 1:100){
    feature_quantiles[j,1,i] = split_q1s[[i]][[3]][j]
    feature_quantiles[j,2,i] = split_medians[[i]][[3]][j]
    feature_quantiles[j,3,i] = split_q3s[[i]][[3]][j]
  }
}

dimnames(feature_quantiles) <- list(1921:2020,
                                    c("Q1", "Median", "Q3"),
                                    feature_names)



plot(1930,1, xlim = c(1921,2020), ylim = range(0,15), xlab="Year", ylab="Range")
for(i in 1:length(feature_names)){
  name = feature_names[i]
  medians = as.numeric(feature_quantiles[,,name][,"Median"])

  avg = mean(medians)
  medians = medians / avg
  
  lines(x=1921:2020, y=medians, col=i)
}

legend("topright", lty=1, col=c(1,2,3,4,5,6,7,8,9,10), legend=feature_names)

Plotting annual median values for the feature for the top 10% of songs by popularity. Standardized such that each feature is on the same scale.

q1s = data.frame(matrix(ncol=3, nrow=0))
medians = data.frame(matrix(ncol=3, nrow=0))
q3s = data.frame(matrix(ncol=3, nrow=0))
feature_names = sort(c("acousticness", "danceability", "instrumentalness", "energy",
                  "duration_ms", "valence", "tempo", "liveness", "loudness", "speechiness"))
for(i in 1921:2020){
  annual_val = subset(df, year==i)
  
  
  top_10 = floor(nrow(annual_val)*0.1)
  annual_val = annual_val[order(annual_val$popularity, decreasing=TRUE),]
  annual_val = annual_val[1:top_10,]

  for(name in feature_names){
    vals = quantile(annual_val[,name])[2:4]
    q1s <- rbind(q1s, c(i, name, vals[1]))
    medians <- rbind(medians, c(i, name, vals[2]))
    q3s <- rbind(q3s, c(i, name, vals[3]))
  }
}

colnames(q1s) <-c("Year", "Feature", "Value")
colnames(medians) <-c("Year", "Feature", "Value")
colnames(q3s) <-c("Year", "Feature", "Value")

split_q1s <- split(q1s, q1s$Feature)
split_medians <- split(medians, medians$Feature)
split_q3s <- split(q3s, q3s$Feature)

feature_quantiles = array(0, c(100,3,length(feature_names)))

for(i in 1:length(feature_names)){
  for(j in 1:100){
    feature_quantiles[j,1,i] = split_q1s[[i]][[3]][j]
    feature_quantiles[j,2,i] = split_medians[[i]][[3]][j]
    feature_quantiles[j,3,i] = split_q3s[[i]][[3]][j]
  }
}

dimnames(feature_quantiles) <- list(1921:2020,
                                    c("Q1", "Median", "Q3"),
                                    feature_names)


plot(1930,1, xlim = c(1921,2020), ylim = range(0,17), xlab="Year", ylab="Range")
for(i in 1:length(feature_names)){
  name = feature_names[i]
  medians = as.numeric(feature_quantiles[,,name][,"Median"])
  avg = mean(medians)
  medians = medians / avg
  
  lines(x=1921:2020, y=medians, col=i)
}

legend("topright", lty=1, col=c(1,2,3,4,5,6,7,8,9,10), legend=feature_names)

Now plotting absolute values for each feature as the median values for overall, as well as just the most popular songs.

MAPE is the difference between the most popular median and the overall median; R is the correlation between overall median and most popular median.

These values will determine whether or not the overall median is a useful predictor of what the median for popular songs will be. This is important because it will be easier to forecast overall values into the future, given they are much smoother, and these forecast overall values can be used to predict the most popular feature value.

# Get the list of features I need to work with

library(Metrics)

q1s = data.frame(matrix(ncol=3, nrow=0))
medians = data.frame(matrix(ncol=3, nrow=0))
best_medians = data.frame(matrix(ncol=3, nrow=0))
q3s = data.frame(matrix(ncol=3, nrow=0))
feature_names = sort(c("acousticness", "danceability", "instrumentalness", "energy",
                  "duration_ms", "valence", "tempo", "liveness", "loudness", "speechiness"))
for(i in 1921:2020){
  annual_val = subset(df, year==i)
  top_10 = floor(nrow(annual_val)*0.1)
  annual_val = annual_val[order(annual_val$popularity, decreasing=TRUE),]
  best_ones = annual_val[1:top_10,]

  for(name in feature_names){
    vals = quantile(annual_val[,name])[2:4]
    best_vals = quantile(best_ones[,name])[2:4]
    best_medians <- rbind(best_medians, c(i, name, best_vals[2]))
    medians <- rbind(medians, c(i, name, vals[2]))
  }
}

colnames(medians) <-c("Year", "Feature", "Value")
colnames(best_medians) <- c("Year", "Feature", "Value")

split_medians <- split(medians, medians$Feature)
split_best_medians <- split(best_medians, best_medians$Feature)

feature_quantiles = matrix(0, ncol=length(feature_names), nrow=100)
best_feature_quantiles = matrix(0, ncol=length(feature_names), nrow=100)

for(i in 1:length(feature_names)){
  for(j in 1:100){
    feature_quantiles[j,i] <- split_medians[[i]][[3]][j]
    best_feature_quantiles[j,i] = split_best_medians[[i]][[3]][j]
  }
}

dimnames(feature_quantiles) <- list(1921:2020,
                                    feature_names)
dimnames(best_feature_quantiles) <- list(1921:2020,
                                    feature_names)




for(i in 1:length(feature_names)){
  name = feature_names[i]
  scores = as.numeric(feature_quantiles[,i])
  best_scores = as.numeric(best_feature_quantiles[,i])
  mape = mape(best_scores, scores)
  mape = round(mape * 100,3)
  r  = cor(scores, best_scores)
  pval = t.test(scores, )
  # Normalize it????? / make it a function of its own mean
  
  smallest = min(min(scores), min(best_scores))
  biggest = max(max(scores), max(best_scores))
  
  plot(x=1921:2020, y=scores, col=1, type="l", ylim=c(smallest, biggest), main=name, sub=paste("MAPE: ", mape, "%;", "R:", r), xlab="Year", ylab="Median Value")
  lines(x=1921:2020, y=best_scores, col=2)
  legend("topright", lty=1, col=c(1,2), legend=c("Overall", "Most Popular"))
  title()
}

Train time series models to forecast the future models.

# Set up the libraries and the training/testing amounts
library("smooth")

Loading required package: greybox
Package "greybox", v0.6.3 loaded.
If you want to know more about the greybox and forecasting, you can visit my website: https://forecasting.svetunkov.ru/


Attaching package: ‘greybox’

The following object is masked from ‘package:lubridate’:

    hm

The following object is masked from ‘package:tidyr’:

    spread

This is package "smooth", v2.6.0

Attaching package: ‘smooth’

The following object is masked from ‘package:TTR’:

    lags

library("forecast")
library("nnfor")
training.percent = 0.95
nTrain = 100*training.percent
nTest = 100*(1-training.percent)

Accousticness Models

acousticness = ts(as.numeric(feature_quantiles[,"acousticness"]) , start=1921)
acousticness.train = subset(acousticness, start=1, end=nTrain)
acousticness.test = subset(acousticness, start = (nTrain+1), end =(nTrain+nTest))

#ses
acousticness.train.ses <- ses(acousticness.train, h=nTest)
acousticness.ses.mape = mape(acousticness.train.ses$mean, acousticness.test)
plot(acousticness.train.ses, main=paste("Acousticness", "SES"), sub=paste("MAPE:", round(acousticness.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(acousticness)

#SARIMA
acousticness.sarima.model <- auto.arima(acousticness.train)
acousticness.sarima <- forecast(acousticness.sarima.model, h=nTest)
acousticness.sarima.mape = mape(acousticness.sarima$mean, acousticness.test)
plot(acousticness.sarima, main=paste("Acousticness", "SARIMA"), sub=paste("MAPE:", round(acousticness.sarima.mape*100,3), "%"), xlab="Year", ylab="Median Value")
lines(acousticness)

#Neural Network
acousticness.train.nn <- elm(acousticness.train)
acousticness.nn.forecast <- forecast(acousticness.train.nn, h=nTest)
acousticness.nn.mape = mape(acousticness.nn.forecast$mean, acousticness.test)
plot(acousticness, main=paste("Acousticness", "NN"), sub=paste("MAPE:", round(acousticness.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(acousticness.nn.forecast$mean, start=1921+nTrain), col=2)

Here we see that the best method to predict future acousticness median values for the overall data set is through SES. However, given that this best MAPE is roughly 37%, this feature cannot be reliably forecast into the future.

Energy models

energy = ts(as.numeric(feature_quantiles[,"energy"]) , start=1921)
energy.train = subset(energy, start=1, end=nTrain)
energy.test = subset(energy, start = (nTrain+1), end =(nTrain+nTest))

#SES
energy.train.ses <- ses(energy.train, h=nTest)
energy.ses.mape = mape(energy.train.ses$mean, energy.test)
plot(energy.train.ses, main=paste("Energy", "SES"), sub=paste("MAPE:", round(energy.ses.mape* 100, 3), "%"), xlab="Year", ylab="Median Value")
lines(energy)

#SARIMA
energy.sarima.model <- auto.arima(energy.train)
energy.sarima <- forecast(energy.sarima.model, h=nTest)
energy.sarima.mape = mape(energy.sarima$mean, energy.test)
plot(energy.sarima, main=paste("Energy", "SARIMA"), sub=paste("MAPE:", round(energy.sarima.mape* 100, 3), "%"), xlab="Year", ylab="Median Value")
lines(energy)

#Neural Network
energy.train.nn <- elm(energy.train)
energy.nn.forecast <- forecast(energy.train.nn, h=nTest)
energy.nn.mape = mape(energy.nn.forecast$mean, energy.test)
plot(energy, main=paste("Energy", "NN"), sub=paste("MAPE:", round(energy.nn.mape* 100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(energy.nn.forecast$mean, start=1921+nTrain), col=2)

Here we see that the best method to predict future energy median values for the overall data set is through SES. Given that this best MAPE is roughly 8%, this feature can be fairly reliably forecast into the future.

Danceability models

danceability = ts(as.numeric(feature_quantiles[,"danceability"]) , start=1921)
danceability.train = subset(danceability, start=1, end=nTrain)
danceability.test = subset(danceability, start = (nTrain+1), end =(nTrain+nTest))

#SES
danceability.train.ses <- ses(danceability.train, h=nTest)
danceability.ses.mape = mape(danceability.train.ses$mean, danceability.test)
plot(danceability.train.ses, main=paste("danceability", "SES"), sub=paste("MAPE:", round(danceability.ses.mape* 100, 3), "%"), xlab="Year", ylab="Median Value")
lines(danceability)

#SARIMA
danceability.sarima.model <- auto.arima(danceability.train)
danceability.sarima <- forecast(danceability.sarima.model, h=nTest)
danceability.sarima.mape = mape(danceability.sarima$mean, danceability.test)
plot(danceability.sarima, main=paste("danceability", "SARIMA"), sub=paste("MAPE:", round(danceability.sarima.mape* 100, 3), "%"), xlab="Year", ylab="Median Value")
lines(danceability)

#Neural Network
danceability.train.nn <- elm(danceability.train)
danceability.nn.forecast <- forecast(danceability.train.nn, h=nTest)
danceability.nn.mape = mape(danceability.nn.forecast$mean, danceability.test)
plot(danceability, main=paste("danceability", "NN"), sub=paste("MAPE:", round(danceability.nn.mape* 100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(danceability.nn.forecast$mean, start=1921+nTrain), col=2)

Here we see that the best method to predict future danceability median values for the overall data set is through NN. However, given that this best MAPE is roughly 10.4%, this feature cannot be reliably forecast into the future.

Duration models

duration = ts(as.numeric(feature_quantiles[,"duration_ms"]) , start=1921)
duration.train = subset(duration, start=1, end=nTrain)
duration.test = subset(duration, start = (nTrain+1), end =(nTrain+nTest))

#SES
duration.train.ses <- ses(duration.train, h=nTest)
duration.ses.mape = mape(duration.train.ses$mean, duration.test)
plot(duration.train.ses, main=paste("Duration", "SES"), sub=paste("MAPE:",round(duration.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(duration)

#SARIMA
duration.sarima.model <- auto.arima(duration.train)
duration.sarima <- forecast(duration.sarima.model, h=nTest)
duration.sarima.mape = mape(duration.sarima$mean, duration.test)
plot(duration.sarima, main=paste("Duration", "SARIMA"), sub=paste("MAPE:", round(duration.sarima.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(duration)

#Neural Network
duration.train.nn <- elm(duration.train)
duration.nn.forecast <- forecast(duration.train.nn, h=nTest)
duration.nn.mape = mape(duration.nn.forecast$mean, duration.test)
plot(duration, main=paste("Duration", "NN"), sub=paste("MAPE:", round(duration.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(duration.nn.forecast$mean, start=1921+nTrain), col=2)

Here we see that the best method to predict future duration median values for the overall data set is through NN. Given that this best MAPE is roughly 8%, this feature can be fairly reliably forecast into the future.

Valence

valence = ts(as.numeric(feature_quantiles[,"valence"]) , start=1921)
valence.train = subset(valence, start=1, end=nTrain)
valence.test = subset(valence, start = (nTrain+1), end =(nTrain+nTest))

#SES
valence.train.ses <- ses(valence.train, h=nTest)
valence.ses.mape = mape(valence.train.ses$mean, valence.test)
plot(valence.train.ses, main=paste("valence", "SES"), sub=paste("MAPE:",round(valence.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(valence)

#SARIMA
valence.sarima.model <- auto.arima(valence.train)
valence.sarima <- forecast(valence.sarima.model, h=nTest)
valence.sarima.mape = mape(valence.sarima$mean, valence.test)
plot(valence.sarima, main=paste("valence", "SARIMA"), sub=paste("MAPE:", round(valence.sarima.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(valence)

#Neural Network
valence.train.nn <- elm(valence.train)
valence.nn.forecast <- forecast(valence.train.nn, h=nTest)
valence.nn.mape = mape(valence.nn.forecast$mean, valence.test)
plot(valence, main=paste("valence", "NN"), sub=paste("MAPE:", round(valence.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(valence.nn.forecast$mean, start=1921+nTrain), col=2)

Tempo model

tempo = ts(as.numeric(feature_quantiles[,"tempo"]) , start=1921)
tempo.train = subset(tempo, start=1, end=nTrain)
tempo.test = subset(tempo, start = (nTrain+1), end =(nTrain+nTest))

#SES
tempo.train.ses <- ses(tempo.train, h=nTest)
tempo.ses.mape = mape(tempo.train.ses$mean, tempo.test)
plot(tempo.train.ses, main=paste("tempo", "SES"), sub=paste("MAPE:",round(tempo.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(tempo)

#SARIMA
tempo.sarima.model <- auto.arima(tempo.train)
tempo.sarima <- forecast(tempo.sarima.model, h=nTest)
tempo.sarima.mape = mape(tempo.sarima$mean, tempo.test)
plot(tempo.sarima, main=paste("tempo", "SARIMA"), sub=paste("MAPE:", round(tempo.sarima.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(tempo)

#Neural Network
tempo.train.nn <- elm(tempo.train)
tempo.nn.forecast <- forecast(tempo.train.nn, h=nTest)
tempo.nn.mape = mape(tempo.nn.forecast$mean, tempo.test)
plot(tempo, main=paste("tempo", "NN"), sub=paste("MAPE:", round(tempo.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(tempo.nn.forecast$mean, start=1921+nTrain), col=2)

Here we see that the best method to predict future tempo median values for the overall data set is through SES. Given that this best MAPE is roughly 1%, this feature can be fairly reliably forecast into the future.

Liveness

liveness = ts(as.numeric(feature_quantiles[,"liveness"]) , start=1921)
liveness.train = subset(liveness, start=1, end=nTrain)
liveness.test = subset(liveness, start = (nTrain+1), end =(nTrain+nTest))

#SES
liveness.train.ses <- ses(liveness.train, h=nTest)
liveness.ses.mape = mape(liveness.train.ses$mean, liveness.test)
plot(liveness.train.ses, main=paste("liveness", "SES"), sub=paste("MAPE:",round(liveness.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(liveness)

#SARIMA
liveness.sarima.model <- auto.arima(liveness.train)
liveness.sarima <- forecast(liveness.sarima.model, h=nTest)
liveness.sarima.mape = mape(liveness.sarima$mean, liveness.test)
plot(liveness.sarima, main=paste("liveness", "SARIMA"), sub=paste("MAPE:", round(liveness.sarima.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(liveness)

#Neural Network
liveness.train.nn <- elm(liveness.train)
liveness.nn.forecast <- forecast(liveness.train.nn, h=nTest)
liveness.nn.mape = mape(liveness.nn.forecast$mean, liveness.test)
plot(liveness, main=paste("liveness", "NN"), sub=paste("MAPE:", round(liveness.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(liveness.nn.forecast$mean, start=1921+nTrain), col=2)

Loudness

loudness = ts(as.numeric(feature_quantiles[,"loudness"]) , start=1921)
loudness.train = subset(loudness, start=1, end=nTrain)
loudness.test = subset(loudness, start = (nTrain+1), end =(nTrain+nTest))

#SES
loudness.train.ses <- ses(loudness.train, h=nTest)
loudness.ses.mape = mape(loudness.train.ses$mean, loudness.test)
plot(loudness.train.ses, main=paste("loudness", "SES"), sub=paste("MAPE:",round(loudness.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(loudness)

#SARIMA
loudness.sarima.model <- auto.arima(loudness.train)
loudness.sarima <- forecast(loudness.sarima.model, h=nTest)
loudness.sarima.mape = mape(loudness.sarima$mean, loudness.test)
plot(loudness.sarima, main=paste("loudness", "SARIMA"), sub=paste("MAPE:", round(loudness.sarima.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(loudness)

#Neural Network
loudness.train.nn <- elm(loudness.train)
loudness.nn.forecast <- forecast(loudness.train.nn, h=nTest)
loudness.nn.mape = mape(loudness.nn.forecast$mean, loudness.test)
plot(loudness, main=paste("loudness", "NN"), sub=paste("MAPE:", round(loudness.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(loudness.nn.forecast$mean, start=1921+nTrain), col=2)

Speechiness

speechiness = ts(as.numeric(feature_quantiles[,"speechiness"]) , start=1921)
speechiness.train = subset(speechiness, start=1, end=nTrain)
speechiness.test = subset(speechiness, start = (nTrain+1), end =(nTrain+nTest))

#SES
speechiness.train.ses <- ses(speechiness.train, h=nTest)
speechiness.ses.mape = mape(speechiness.train.ses$mean, speechiness.test)
plot(speechiness.train.ses, main=paste("speechiness", "SES"), sub=paste("MAPE:",round(speechiness.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(speechiness)

#SARIMA
speechiness.sarima.model <- auto.arima(speechiness.train)
speechiness.sarima <- forecast(speechiness.sarima.model, h=nTest)
speechiness.sarima.mape = mape(speechiness.sarima$mean, speechiness.test)
plot(speechiness.sarima, main=paste("speechiness", "SARIMA"), sub=paste("MAPE:", round(speechiness.sarima.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(speechiness)

#Neural Network
speechiness.train.nn <- elm(speechiness.train)
speechiness.nn.forecast <- forecast(speechiness.train.nn, h=nTest)
speechiness.nn.mape = mape(speechiness.nn.forecast$mean, speechiness.test)
plot(speechiness, main=paste("speechiness", "NN"), sub=paste("MAPE:", round(speechiness.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(speechiness.nn.forecast$mean, start=1921+nTrain), col=2)

Linear Model for Overall Feature Median–>Most Popular Songs’ Feature Median

We will only consider models that had a low enough MAPEs (10% chosen as a threshold) in their best timeseries forecasting model.

Here MAPE is the difference between the actual median for most popular songs and the predicted median.

#Energy
best_energy = ts(as.numeric(best_feature_quantiles[,"energy"]), start=1921, end=2020) 
best_energy.train = subset(best_energy, start = 1, end = nTrain)
best_energy.test = subset(best_energy, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_energy.train, overall = energy.train)

best = best_energy.train
overall = energy.train

energy.lm <- lm(best ~ overall)
test_df <- data.frame(overall = energy.test)
best_energy.test.predict <- predict(energy.lm, newdata=test_df, interval="prediction")
best_energy_test_predict.mape <- mape(best_energy.test[1:length(best_energy.test)], best_energy.test.predict[,1])
plot(best_energy,xlab="Year", ylab="Energy", main=paste("Energy Popularity Prediction"), sub = paste("MAPE: ", round(best_energy_test_predict.mape*100, 3), "%"))
lines(ts(best_energy.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Duration
best_duration = ts(as.numeric(best_feature_quantiles[,"duration_ms"]), start=1921, end=2020) 
best_duration.train = subset(best_duration, start = 1, end = nTrain)
best_duration.test = subset(best_duration, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_duration.train, overall = duration.train)

best = best_duration.train
overall = duration.train

duration.lm <- lm(best ~ overall)
test_df <- data.frame(overall = duration.test)
best_duration.test.predict <- predict(duration.lm, newdata=test_df, interval="prediction")
best_duration_test.predict.mape <- mape(best_duration.test[1:length(best_duration.test)], best_duration.test.predict[,1])
plot(best_duration,xlab="Year", ylab="duration", main=paste("duration Popularity Prediction"), sub = paste("MAPE: ", round(best_duration_test.predict.mape*100, 3), "%"))
lines(ts(best_duration.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Accousticness
best_acousticness = ts(as.numeric(best_feature_quantiles[,"acousticness"]), start=1921, end=2020) 
best_acousticness.train = subset(best_acousticness, start = 1, end = nTrain)
best_acousticness.test = subset(best_acousticness, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_acousticness.train, overall = acousticness.train)

best = best_acousticness.train
overall = acousticness.train

acousticness.lm <- lm(best ~ overall)
test_df <- data.frame(overall = acousticness.test)
best_acousticness.test.predict <- predict(acousticness.lm, newdata=test_df, interval="prediction")
best_acousticness_test.predict.mape <- mape(best_acousticness.test[1:length(best_acousticness.test)], best_acousticness.test.predict[,1])
plot(best_acousticness,xlab="Year", ylab="acousticness", main=paste("acousticness Popularity Prediction"), sub = paste("MAPE: ", round(best_acousticness_test.predict.mape*100, 3), "%"))
lines(ts(best_acousticness.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Valence
best_valence = ts(as.numeric(best_feature_quantiles[,"valence"]), start=1921, end=2020) 
best_valence.train = subset(best_valence, start = 1, end = nTrain)
best_valence.test = subset(best_valence, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_valence.train, overall = valence.train)

best = best_valence.train
overall = valence.train

valence.lm <- lm(best ~ overall)
test_df <- data.frame(overall = valence.test)
best_valence.test.predict <- predict(valence.lm, newdata=test_df, interval="prediction")
best_valence_test.predict.mape <- mape(best_valence.test[1:length(best_valence.test)], best_valence.test.predict[,1])
plot(best_valence,xlab="Year", ylab="Valence", main=paste("Valence Popularity Prediction"), sub = paste("MAPE: ", round(best_valence_test.predict.mape*100, 3), "%"))
lines(ts(best_valence.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Tempo
best_tempo = ts(as.numeric(best_feature_quantiles[,"tempo"]), start=1921, end=2020) 
best_tempo.train = subset(best_tempo, start = 1, end = nTrain)
best_tempo.test = subset(best_tempo, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_tempo.train, overall = tempo.train)

best = best_tempo.train
overall = tempo.train

tempo.lm <- lm(best ~ overall)
test_df <- data.frame(overall = tempo.test)
best_tempo.test.predict <- predict(tempo.lm, newdata=test_df, interval="prediction")
best_tempo_test.predict.mape <- mape(best_tempo.test[1:length(best_tempo.test)], best_tempo.test.predict[,1])
plot(best_tempo,xlab="Year", ylab="tempo", main=paste("tempo Popularity Prediction"), sub = paste("MAPE: ", round(best_tempo_test.predict.mape*100, 3), "%"))
lines(ts(best_tempo.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Liveness
best_liveness = ts(as.numeric(best_feature_quantiles[,"liveness"]), start=1921, end=2020) 
best_liveness.train = subset(best_liveness, start = 1, end = nTrain)
best_liveness.test = subset(best_liveness, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_liveness.train, overall = liveness.train)

best = best_liveness.train
overall = liveness.train

liveness.lm <- lm(best ~ overall)
test_df <- data.frame(overall = liveness.test)
best_liveness.test.predict <- predict(liveness.lm, newdata=test_df, interval="prediction")
best_liveness_test.predict.mape <- mape(best_liveness.test[1:length(best_liveness.test)], best_liveness.test.predict[,1])
plot(best_liveness,xlab="Year", ylab="liveness", main=paste("liveness Popularity Prediction"), sub = paste("MAPE: ", round(best_liveness_test.predict.mape*100, 3), "%"))
lines(ts(best_liveness.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Loudness
best_loudness = ts(as.numeric(best_feature_quantiles[,"loudness"]), start=1921, end=2020) 
best_loudness.train = subset(best_loudness, start = 1, end = nTrain)
best_loudness.test = subset(best_loudness, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_loudness.train, overall = loudness.train)

best = best_loudness.train
overall = loudness.train

loudness.lm <- lm(best ~ overall)
test_df <- data.frame(overall = loudness.test)
best_loudness.test.predict <- predict(loudness.lm, newdata=test_df, interval="prediction")
best_loudness_test.predict.mape <- mape(best_loudness.test[1:length(best_loudness.test)], best_loudness.test.predict[,1])
plot(best_loudness,xlab="Year", ylab="loudness", main=paste("loudness Popularity Prediction"), sub = paste("MAPE: ", round(best_loudness_test.predict.mape*100, 3), "%"))
lines(ts(best_loudness.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

Predict features for the overall market trends for the next five years.

Only including features whose trained linear models had a MAPE below 10%, as these are the features whose most popular medians will be most reliably predicted by their forecast overall values.

Using the best timeseries method seen for each feature.

Energy: SES Model

energy.ses <- ses(energy, h=nTest)
plot(energy.ses, main=paste("Energy", "SES"), xlab="Year", ylab="Median Value")
lines(energy)

Duration: NN Model

duration.nn <- elm(duration)
duration.nn.predict <- forecast(duration.nn, h=nTest)
plot(duration.nn.predict, main=paste("Duration", "NN"), xlab="Year", ylab="Median Value")

Tempo: SES Model

tempo.ses <- ses(tempo, h=nTest)
plot(tempo.ses, main=paste("Tempo", "SES"), xlab="Year", ylab="Median Value")
lines(tempo)

Valence: SARIMA Model

valence.sarima <- auto.arima(valence)
valence.sarima.predict <- forecast(valence.sarima, h=nTest)
plot(valence.sarima.predict, main=paste("Valence", "SARIMA"), xlab="Year", ylab="Median Value")
lines(valence)

Liveness: SARIMA Model

liveness.sarima <- auto.arima(liveness)
liveness.sarima.predict <- forecast(liveness.sarima, h=nTest)
plot(liveness.sarima.predict, main=paste("liveness", "SARIMA"), xlab="Year", ylab="Median Value")
lines(liveness)

Loudness: SES Model

loudness.ses <- ses(loudness, h=nTest)
plot(loudness.ses, main=paste("loudness", "SES"), xlab="Year", ylab="Median Value")
lines(loudness)

Now predicting the median values for the most popular songs for this set of features.

The forecast overall median values and linear model predicting the most popular median values for each feature given overall median value will be used to predict the future median values for popular songs for each feature.

Energy

best = best_energy
overall = energy


energy.lm <- lm(best ~ overall)
predict_df <- data.frame(overall = energy.ses$mean)
energy.predict <- predict(energy.lm, newdata=predict_df, interval="prediction")

smallest = min(min(best), min(overall))
biggest = max(max(best), max(overall))
plot(energy, xlab="Year", ylab="Median Value", ylim=c(smallest, biggest), main="Most Popular Song's Energy Feature Prediction", sub=paste("Projected median feature value for most popular songs:", paste(round(energy.ses$mean, 3), collapse=",") ))
lines(best_energy, col = 2)

lines(ts(energy.predict[,1], start=2020+1, end=2020+nTest), col = 3)
legend("bottomright", lty=1, col=c(1,2,3), legend=c("Overall", "Most Popular", "Most Popular Projection"))

Duration

best = best_duration
overall = duration

duration.lm <- lm(best ~ overall)
predict_df <- data.frame(overall = duration.nn.predict$mean)
duration.predict <- predict(duration.lm, newdata=predict_df, interval="prediction")

smallest = min(min(best), min(overall))
biggest = max(max(best), max(overall))
plot(duration, xlab="Year", ylab="Median Value", ylim=c(smallest, biggest), main="Most Popular Song's Duration Feature Prediction", sub=paste("Projected median feature value for most popular songs:", paste(round(duration.nn.predict$mean), collapse=",") ))
lines(best_duration, col = 2)

lines(ts(duration.predict[,1], start=2020+1, end=2020+nTest), col = 3)
legend("topright", lty=1, col=c(1,2,3), legend=c("Overall", "Most Popular", "Most Popular Projection"))

Tempo

best = best_tempo
overall = tempo

tempo.lm <- lm(best ~ overall)
predict_df <- data.frame(overall = tempo.ses$mean)
tempo.predict <- predict(tempo.lm, newdata=predict_df, interval="prediction")


smallest = min(min(best), min(overall))
biggest = max(max(best), max(overall))
plot(tempo, xlab="Year", ylab="Median Value", ylim=c(smallest, biggest), main="Most Popular Song's Tempo Feature Prediction", sub=paste("Projected median feature value for most popular songs:", paste(round(tempo.ses$mean), collapse=",") ))
lines(best_tempo, col = 2)

lines(ts(tempo.predict[,1], start=2020+1, end=2020+nTest), col = 3)
legend("bottomright", lty=1, col=c(1,2,3), legend=c("Overall", "Most Popular", "Most Popular Projection"))

Valence

best = best_valence
overall = valence

valence.lm <- lm(best ~ overall)
predict_df <- data.frame(overall = valence.sarima.predict$mean)
valence.predict <- predict(valence.lm, newdata=predict_df, interval="prediction")


smallest = min(min(best), min(overall))
biggest = max(max(best), max(overall))
plot(valence, xlab="Year", ylab="Median Value", ylim=c(smallest, biggest), main="Most Popular Song's valence Feature Prediction", sub=paste("Projected median feature value for most popular songs:", paste(round(valence.sarima.predict$mean,3), collapse=",") ))
lines(best_valence, col = 2)

lines(ts(valence.predict[,1], start=2020+1, end=2020+nTest), col = 3)
legend("bottomright", lty=1, col=c(1,2,3), legend=c("Overall", "Most Popular", "Most Popular Projection"))

Liveness

best = best_liveness
overall = liveness

liveness.lm <- lm(best ~ overall)
predict_df <- data.frame(overall = liveness.sarima.predict$mean)
liveness.predict <- predict(liveness.lm, newdata=predict_df, interval="prediction")

smallest = min(min(best), min(overall))
biggest = max(max(best), max(overall))
plot(liveness, xlab="Year", ylab="Median Value", ylim=c(smallest, biggest), main="Most Popular Song's Liveness Feature Prediction", sub=paste("Projected median feature value for most popular songs:", paste(round(liveness.sarima.predict$mean,3), collapse=",") ))
lines(best_liveness, col = 2)

lines(ts(liveness.predict[,1], start=2020+1, end=2020+nTest), col = 3)
legend("topright", lty=1, col=c(1,2,3), legend=c("Overall", "Most Popular", "Most Popular Projection"))

Loudness

best = best_loudness
overall = loudness

loudness.lm <- lm(best ~ overall)
predict_df <- data.frame(overall = loudness.ses$mean)
loudness.predict <- predict(loudness.lm, newdata=predict_df, interval="prediction")

smallest = min(min(best), min(overall))
biggest = max(max(best), max(overall))
plot(loudness, xlab="Year", ylab="Median Value", ylim=c(smallest, biggest), main="Most Popular Song's Loudness Feature Prediction", sub=paste("Projected median feature value for most popular songs:", paste(round(loudness.ses$mean,1), collapse=",") ))
lines(best_loudness, col = 2)

lines(ts(loudness.predict[,1], start=2020+1, end=2020+nTest), col = 3)
legend("bottomright", lty=1, col=c(1,2,3), legend=c("Overall", "Most Popular", "Most Popular Projection"))

Implications:

Using these values above for the predicted median value for the most popular songs, for each feature, music curators and record labels can predict which “sounds” will be most popular in the future. This information can then be used for future decision making with respect to which artists to sign and which sounds to push.